% The code is based on the following paper:
%   Chen, T. Y., Lin, Y. L., and Tzeng, L. Y., forthcoming.
%   Estimating probability weighting functions through option pricing bounds. Review of Asset Pricing Studies.
%
% Copyright: Tzu-Ying Chen, Yo-Lan Lin, Larry Y. Tzeng
% Date: January 30, 2024

clear

%% Setting                
BV_Target = {'UB'};                                                        
             
%% Optimal Risk Preferences under Rank-Dependent Expected Utility
FileName = ['Summary_OPT_AllCR_AllPeriod_OP_UB'];
load(FileName, ...
     'Date_Monthly', ...
     'Index_CR', 'Index_OPT', 'OPT_TAB');
clear FileName 

OPT = nan(size(Date_Monthly, 1), length(Index_OPT));                       
for d = 1:size(Date_Monthly, 1)
    OPT(d, :) = OPT_TAB(Index_CR, Index_OPT, d);                                 
end 
clear Index_CR Index_OPT OPT_TAB
clear d

%% Mean and Variance of the Return Distribution
AllMoment_Text = {'$\mu_{t} = E^{v}_{H} \left[ S_{T} / S_{t} \mid S_{t} \right]$'; ...
                  '$\sigma_{t}^{2} = E^{v}_{H} \left[ \left( S_{T} / S_{t} - \mu_{t} \right)^{2} \mid S_{t} \right]$'};
      
AllMoments_EU = nan(size(Date_Monthly, 1), size(AllMoment_Text, 1));       
AllMoments_RDEU = nan(size(Date_Monthly, 1), size(AllMoment_Text, 1));     
for d = 1:size(Date_Monthly, 1)
    Target_Date = Date_Monthly(d, 1);
    Target_TTM = Date_Monthly(d, 3);
    
    % Load Data 
    FileName = ['Summary_OP_Grid_' num2str(Target_Date)];
    load(FileName, ...
         'State_Monthly', 'State_PROB');
    clear FileName  

    % Expected Utility
    M1 = sum(State_PROB .* power(State_Monthly, 1));
    AllMoments_EU(d, 1) = M1;
    AllMoments_EU(d, 2) = sum(State_PROB .* power(State_Monthly - M1, 2));
    clear M1
    
    % Rank-Dependent Expected Utility
    CDF = cumsum(State_PROB);
    CDF(end) = 1;                                                           
    CDF = P1998(CDF, OPT(d, 1), OPT(d, 2));                                
    State_PROB = [CDF(1) diff(CDF)];                                         
    clear CDF

    M1 = sum(State_PROB .* power(State_Monthly, 1));
    AllMoments_RDEU(d, 1) = M1;
    AllMoments_RDEU(d, 2) = sum(State_PROB .* power(State_Monthly - M1, 2));
    clear Target_Date Target_TTM State_Monthly State_PROB M1    
end
clear d

%% Plot Figure: Mean and Variance of the Return Distribution
AllDate = datenum(num2str(Date_Monthly(:, 1)), 'yyyymmdd');               

h_Moments = figure;
for i = 1:size(AllMoment_Text, 1)
    subplot(2, 1, i)
    plot(AllDate, AllMoments_RDEU(:, i), ...
         'LineStyle', '-', ...
         'LineWidth', 2, ...
         'Color', 'r')     
    hold on      
 
    plot(AllDate, AllMoments_EU(:, i), ...
         'LineStyle', ':', ...
         'LineWidth', 2, ...
         'Color', 'b')     
    grid on         
    
    h = legend(['($\alpha$ , $\beta$) = ($\alpha_{t}^{\ast}$ , $\beta_{t}^{\ast}$)'], ...
               ['($\alpha$ , $\beta$) = (' num2str(1) ' , ' num2str(1) ')']);
    set(h, 'FontSize', 15, ...
           'FontName', 'Times New Roman', ...
           'FontWeight', 'Bold', ...
           'Location', 'NorthEast', ...
           'Interpreter', 'Latex', ...
           'Box', 'Off')
    clear h

    h1 = xlabel('Time');  
    h2 = ylabel('Value');
    h3 = title(AllMoment_Text(i, :));
    set([h1 h2 h3], 'FontSize', 15, ...
                    'FontName', 'Times New Roman', ...
                    'FontWeight', 'Bold', ...
                    'Interpreter', 'Latex')
    clear h1 h2 h3       
    
    set(gca, 'XTick', [AllDate(1); AllDate(52:52:end)], ...
             'FontSize', 15, ...
             'FontName', 'Times New Roman')
    datetick('x', 'mmmyyyy', 'KeepTicks')   
      
    if i==1
        set(gca, 'YLim', [0.3 1.4], ...
                 'YTick', 0:0.1:1.5, ...
                 'FontSize', 15, ...
                 'FontName', 'Times New Roman')          
    else
    end
    
    if i==2
        set(gca, 'YLim', [0 0.13], ...
                 'YTick', 0:0.01:1.5, ...
                 'FontSize', 15, ...
                 'FontName', 'Times New Roman')          
    else
    end    
    
    set(gca, 'Layer', 'Top')
    set(gca, 'LooseInset', get(gca, 'TightInset'))   
    set(gcf, 'Position', get(0, 'ScreenSize'))
end
clear i

X = 29.7;                                                                  
Y = 21.0;                                                                  
XMargin = 1;                                                               
YMargin = 1;                                                               
XSize = X - 1 * XMargin;                                                   
YSize = Y - 1 * YMargin;                                                   
set(h_Moments, 'PaperUnits', 'Centimeters')
set(h_Moments, 'PaperSize', [X Y])
set(h_Moments, 'PaperPosition', [XMargin YMargin XSize YSize])
set(h_Moments, 'PaperOrientation', 'Portrait')
clear AllMoment_Text AllMoments_RDEU AllMoments_EU AllDate h_Moments X Y XMargin YMargin XSize YSize
clear Date_Monthly

%% Output
FileName = ['Figure_DISTMoments_Monthly_UB'];
saveas(gcf, [FileName '.eps'], 'epsc')
print(FileName, '-dpdf', '-fillpage')
winopen([FileName '.pdf'])
clear FileName
